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The use of linear interpolation and sim- 
plified iteration procedures are shown to 
introduce inaccuracies to the rectangu- 
lar grid method of characteristics, par- 
ticularly when applied to subcritical 
flows. Comparisons of experimental and 
computational results are presented illus- 
trating the use of Everett and Newton- 
Gregory interpolation, in addition to a 
more complex iteration procedure, to 
substantially improve the method's abil- 
ity to maintain both steady uniform 
flows under subcritical conditions, and 
retain wave steepness during propaga- 



tion along the drainage pipe. The results 
presented will be directly applicable to 
the building drainage network model 
previously developed at Brunei Univer- 
sity with the support of NBS CBT grant 
aid. 
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1. Introduction 



A computer model that will simulate the propa- 
gation of steep fronted waves along a drainage pipe 
is being developed from recent research by Bridge 
[l] 1 and Swaffield et al. [2]. The rectangular grid 
(R.G.) method of characteristics, employed to 
solve numerically the appropriate unsteady flow 
equations, was used to model the attenuation of the 
flow profiles in the pipe. Several papers describe 
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the derivation and applications of the numerical so- 
lutions of the governing equations for transient par- 
tially filled pipe flow [1,2,3]. This paper presents 
results of an investigation into the computational 
problems identified by Swaffield and Bridge. 
One example is the collapse of the water surface 
profile over time during the passage of steady 
flows under subcritical conditions. It was felt that 
these inaccuracies were probably caused by the lin- 
ear interpolation and simplified iteration proce- 
dures being used within the method of 
characteristics. 

Numerous computer runs were made to deter- 
mine if the accuracy of the method could be im- 
proved by employing a nonlinear interpolation 
procedure and by using a more complex iteration 
procedure. The computed results were compared 
to data from experiments carried out on the Brunei 
test rig. 
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NOTATION 


A 


= water depth cross-sectional area (m 2 ) 


c 


= local wavespeed (m/s) 


D 


= pipe diameter (m) 


/» 


=/(*») 


8 


= acceleration due to gravity (m/s 1 ) 


h 


= constant step length (m) 


So 


=pipe slope 


s t 


= slope of energy grade line 


T 


= water surface width (m) 


t 


= time (s) 


TFAC = time factor 


V 


= local flow average velocity (m/s) 


X 


= distance from upstream boundary (m) 


y 


= local flow depth (m) 


B'/» 


= i th central difference at x =x„. 



2. The Rectangular Grid Method 
of Characteristics 

The St. Venant equations are a pair of quasi-lin- 
ear hyperbolic partial differential equations. The 
equation of continuity is expressed 



v£ + a£ + t2-. 



and the dynamic equation, 



g(S t -S a )+g^ + V~ +~=0. 



dx 



(1) 



(2) 



The St. Venant equations can be transformed into 
their characteristics in a number of ways; Fox [3] 
gives one method. The equations describing the 
positive and negative characteristics are: 



£+i*+«»-«-o 



dx 
dr 



= V+c 



dt c dt +g ^ ^ U 



dx 
dt 



= V~c 



(4) 
(5) 
(6) 



where c~g VA/T. The positive characteristic is 
described by eqs (3 and 4) and the negative charac- 
teristic by (5 and 6). 

These equations are then numerically integrated. 
An approximation is introduced in the integration 
of some terms as the variation of V and y is not 
known as a function of time. Referring to figure 1, 
if the velocity and depth are known at points R and 
S at time level t - At then a first-order approxima- 
tion leads to the following set of equations which 
may be written in terms of the unknown depth and 
velocity at point P at time level t, 

V P ~V R +f- (y e -y R )+g (S R -S a ) At =0 (7) 



x P -x SL =(V R +c B )At 

V P ~ V s +& fyp-y s )+g (S s -So) At =0 

x P ~x s =(V s -c$)At . 



(8) 
(9) 

(10) 

These equations are paired such that eq (7) is true 
only if eq (8) is satisfied and similarly for eqs (9 and 
10). In most cases a first order integration is satis- 
factory, however attention must be given to the 
choice of the time step to ensure a stable solution, 
as mentioned in [3]. Generally, applying the 
Courant condition, 



Af = 



Ax 



(K+c)„ 



(11) 



(3) 



ensures that the points R and S fall within ± Ax of 
point P. Since the velocity and celerity both vary 
throughout the duration of the analysis the maxi- 
mum value of each is found at every time step. 

The algebraic solution of eqs (7-10) reduces to a 
problem of determining the values of jcr, Yh, Vr, Xs, 
y s and V s since * p , f p , t R and t s are known values. 
Referring to figure 1, x R and jc s are determined by 
using eqs (4 and 6) respectively to estimate the t-At 
line intersections of the positive and negative char- 
acteristtics passing through point P. The quantities 
Vr, ys, V R and V s are then determined by interpola- 
tion between grid points A and C and C and B 
respectively for this example of subcritical flow. 
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S B 
Fr - Froude number 

If Fr< 1 then flow is subcritical and S lies between C and B. 
If Fr?l then flow is supercritical and S' lies between ft and C, 



Figure 1- Specified time interval 
technique. 



Finally the solution to point P is obtained through 
the simultaneous solution of eqs (7 through 10). 

3. Inaccuracies in the Method 
of Characteristics 

A summary of the initial and boundary condi- 
tions required to apply the R.G. method of charac- 
teristics to drainage flow is presented in [1], The 
method described above was used in a model simu- 
lating the movement of a hydraulic jump in a 
drainage pipe. Like Streeter and Wylie [4], the po- 
sitions of x R and x. were determined by assuming 
that the slope of the characteristic lines at P were 
equal to those at C. The depth and velocity at 
points R and S were determined by linearly inter- 
polating between grid points A and C and C and B 
(subcritical flow). However, it was found that two 
problems arose. The first being that the water sur- 
face profile downstream of the jump collapsed 
over time, and the second that the hydraulic jump 
acquired a positive velocity during the passage of 
steady flows. Harris [5] mentions that the Streeter- 
Wylie approximations are reasonable provided that 
the values of v and V at A, B and C do not differ 
greatly. The values of depth and velocity across a 
hydraulic jump or a steep-fronted wave do differ 
greatly; therefore, these approximations are not 
valid for the case being considered and a more ex- 
act method must be employed. 



4. Improvement of the R.G. Method 
of Characteristics 

4.1 Integration 

There are several ways by which the R.G. 
method of characteristics can be made more accu- 
rate. One method [6] is to apply the trapezoidal 
rule when integrating the characteristic equations. 
However, the resulting equations must be solved 
by an implicit method. The implicit scheme is fairly 
complicated for open channel flow problems and is 
seldom used as it requires an iteration procedure. 
In addition, the benefits derived from the more ac- 
curate integration of the characteristic equations 
are often lost in the interpolation methods required 
to evaluate conditions at R and S. 

4.2 Iteration 

Accuracy can also be increased by using an itera- 
tive method to determine the positions of ;t R and x$ 
from the slopes of the characteristics. A simplified 
iterative method similar to one developed by Lister 
[7] has been used to improve the computer model. 
The position of jc r is initially determined by assum- 
ing that the slope of the characteristic line at P is 
equal to that at C;y R and V R are then found from a 
nonlinear interpolation of values at surrounding 
gridpoints. (The exact interpolation method will be 
discussed below.) x R is then recalculated assuming 
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that the slope of the characteristic line from P to R 
is equal to that of the previous value of x R : 



i,-i R M=(H+^d( 



(12) 



This process is repeated until the difference be- 
tween successive iterations of x R is less than 0.001 
meters. The same process is used to find jc s . 

43 Interpolation 

The inaccuracies caused by the effect of linearly 
interpolating between grid points to find the values 
of the depth and velocity at points R and S at each 
time step are illustrated in figure 2. The figure 
shows a section of backwater profile at the end of a 
pipe in which the flow is sub critical. Alternatively, 
this can be considered as the backwater profile be- 
tween a hydraulic jump and a junction. For steady 
uniform flow the backwater profile should remain 
constant. However, linear interpolation between 
points A and C and points B and C produce values 
for depth at points R and S which are lower than 
the actual values. This leads to an underestimation 
of the depth at P at the next time step, which in 
turn leads to an overestimation of the velocity at P. 
These interpolation errors build up as the calcula- 
tions continue in time causing the backwater pro- 
file to collapse and the discharge to increase. These 
effects have also been noted by Fox (see [1]). 

Quadratic interpolation [7-9] and Langrangian 
[10-12] interpolation procedures have been used 



in the past to increase the accuracy of the R.G. 
method of characteristics. Jolly and Vevjevich [10] 
found that the Langrangian interpolation yielded 
more accurate results than linear interpolation. It 
was found that accuracy was not improved by us- 
ing a more sophisticated scheme than third-order 
Langrangian interpolation equations combined 
with an iterative procedure to find the positions of 
RandS. 

Although the Larjgrange method minimizes the 
amount of calculations it has the disadvantage that 
it gives no indication of the accuracy available. As 
a result, difference methods tend to be used in ex- 
ploratory work and the Langrange methods in 
well-undei stood routine work [13]. In general, the 
error in using a difference method is of the order of 
the first term to be omitted. Formulae which in- 
volve central differences give the best accuracy for 
a given amount of computation. They also have the 
philosophical merit of giving equal weight to the 
tabulated data on each side of the relevant data 
point. The practical interpolation formulae are 
those of Everett and BesseL If third differences are 
negligible for the accuracy required, then Bessel's 
formula is simpler and hence better. If third differ- 
ences are greater than 60 units in the last decimal 
they cannot be neglected. Thus Everett's formula 
is best for our purposes [14]. Everett's formula in- 
volves only even-order central differences: 




LINEAR INTERPOLATION 



Figure 2-Interpolation error on 
the backwater profile. 



STEADY STATE FLOW 



152 



Journal of Research of the National Bureau of Standards 



Y(x)=Eof +E 2 h 2 f ±E$%+ ... 

+F a f 1 +F 2 h 2 fi+FJ> A fi+ ... 
where 

£o=1-h 

_. -u(l~u)(2-u) 

E 2 ~ 3i 

_ -(-l-u)(u)g~u)(2-u)(3-u) 
A 5! 

' 2 3! 



(13) 

(14) 
(15) 

(16) 
(17) 
(18) 

(19) 

(20) 

Near the ends of the pipe where central differences 
cannot be obtained die Newton-Gregory forward 
backward difference formulae are used. The New- 
ton-Gregory forward difference formula is ex- 
pressed by: 

P(x)=-p(x +hu)=f +uAf a +l [u(u -1)A%] 



„ u(u 2 ~l)(u 2 -4) 
4 5! 



and 



„ == k_£fll j 0<u<1 



+ 



|[»(«-l)(«-2)A 3 /o]- 



...+£[»(»-!). ..(»-n+l)A?a.] (21) 
and the backward difference formula by: 
P(*) = P(* +/i«)=/o+HV/ +! [»(« +l)V 2 / ] 
+ |[«(« + 1)(h + 2)V 3 / ]+... 
. . . +i [«(« + 1) . . . (h +« - 1)V% ]. (22) 

5. Comparison of the Results 

An example of the effects of linear interpolation 
can be seen in figure 3 which shows subcritical par- 
tially filled flow with a wave entering the pipe at a 
time of 0.23 seconds. As time progressed the back- 
water profile collapsed. A time factor, TFAC, of 
5.0 was used in this example. The use of a time 
factor causes the time step eq (11) to become 
smaller: 



Af = 



Ax 



(V^+c^TFAC 



(23) 



It is important to determine a value of TFAC that 
results in a converged solution, whilst not requiring 
excessive computation time. If TFAC is chosen suf- 
ficiently large then a converged solution will al- 
ways be achieved, provided that the results are not 
dominated by rounding and truncation errors. Gen- 
erally, the smallest TFAC that will allow stability is 



depth 

(mm) 



~i 1 — r T"^t — i r 



i ! i r 



-| 1 ) r 



LINEAR INTERPOLATION 
L - 15m 
TFAC - 5.0 
E, - V3QQ 
□ ** O.la 




\t ^ 0,23 sees 



Figure 3. Linear interpolation used when a wave enters a pips in 
which thexe is subcritical flow 
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Figure 3-Linear interpolation used 
when a wave enters a pipe in 
which there is subcritical flow. 
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used. The time factors used in these results were 
determined by trial and error. 

Figure 4 shows the effects of using Everett and 
Newton-Gregory interpolation on the same flow 
conditions of the previous example. A TFAC of 1.0 
is the smallest TFAC that would allow stability, 
from eq (11). It can be seen that the backwater 
profile remained unchanged with increasing time. 
On comparing the two figures it can be seen that 
the linear method tended to flatten out the wave 
whereas Everett's method retained the wave's 



steepness. Everett's method also remained stable 
with a TFAC of 1.0. When linear interpolation was 
used with a TFAC of 1.0 a discontinuity developed 
on the wave front at a time of 11.31 seconds as 
shown in figure 5. 

Figure 6 shows the experimentally measured 
depth profile at a distance of 8.0 meters from the 
entrance of the pipe against time. It can be seen 
that in this case both Everett's method and the lin- 
ear interpolation method used with a TFAC of 1.0 
retained the wave's steepness and compare most 
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Figure 4, Everett and Newton-Gregory interpolation used on 
some flow conditions as Figure 3 
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Figure 4-Everett and Newton- 
Gregory interpolation used on 
some flow conditions as figure 
3. 
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Figure 5. Discontinuity in wave front when using linear interpolation 
with small TFAC 



I I I I I I 1 I ! 1 1 1 1 L. 
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Figure 6-Comparison of experimental and computed results used to predict wave attenuation. 



favorably with the experimental data from the 
Brunei test rig. 

An increasing backwater profile is created 
where a drainage pipe joins a junction. Figure 7 
shows the propagation of a wave entering such a 
pipe in a case where Everett and Newton-Gregory 



interpolation with a TFAC of 3.0 was used. It can 
be seen that the backwater profile remained steady 
and the steepness of the wave was retained. When 
linear interpolation was used with the same TFAC 
the backwater profile increased as time progressed. 
Thus Everett's method is clearly an improvement 
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Figure7 . Everett and Hew tan-Gregory interpolation used on. an Increasing 
backwater profile 
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Figure 7-Everett and Newton- 
Gregory interpolation used on 
an increasing backwater profile. 
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to linear interpolation as it retains the backwater 
profile and wave steepness when run with a small 
TFAC. 



6. Conclusions 

A computer model, based on the rectangular 
grid method of characteristics that will simulate 
the propagation of steep fronted waves along a 
drainage pipe is being developed from recent re- 
search by Bridge and Swaffield [1,2]. Experimental 
tests carried out on the Brunei test rig as well as 
numerous computer runs have shown that the ac- 
curacy of the model has been considerably im- 
proved by including an interation procedure to find 
the positions of the space time coordinates R and S 
and also by using Everett's interpolation to deter- 
mine the height and velocity at points R and S. The 
improvements made have solved the problem of 
the decreasing backwater profile during steady uni- 
form flow noted by Bridge. The improvements also 
retain the steepness of waves shown in test data as 
they travel along the branch pipe. 
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